Tracking shocked dust: state estimation for a complex plasma during a shock wave 
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We consider a two-dimensional complex (dusty) plasma crystal excited by an electrostatically-induced shock 
wave. Dust particle kinematics in such a system are usually determined using particle tracking velocimetry. 
In this work we present a particle tracking algorithm which determines the dust particle kinematics with 
significantly higher accuracy than particle tracking velocimetry. The algorithm uses multiple extended Kalman 
filters to estimate the particle states and an interacting multiple model to assign probabilities to the different 
filters. This enables the determination of relevant physical properties of the dust, such as kinetic energy and 
kinetic temperature, with high precision. We use a Hugoniot shock-jump relation to calculate a pressure- 
volume diagram from the shocked dust kinematics. Calculation of the full pressure-volume diagram was 
possible with our tracking algorithm, but not with particle tracking velocimetry. 
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I. INTRODUCTION 

A complex plasma consists of mesoscopic particles of 
'dust' suspended within an ionized gas — a low-density 
plasma consisting of neutral atoms, ions and electrons. 1 ' 2 
The dust particles repel each other due to acquiring a 
net negative charge from collisions with electrons in the 
plasma (and less frequently, ions). Confined dust can 
form ordered, crystal-like structures. 3 ' 4 In dusty plasma 
monolayer experiments, single particles of dust are re- 
solvable using a digital camera and illumination by a 
two-dimensional (2D) laser sheet, making dusty plasma 
monolayers a useful macromolecular-like system for ex- 
ploring the microscopic kinematics of condensed-matter 
systems. The setup shown in Fig. 5 in Sec. IV can be used 
in experiments to generate and observe a range of phe- 
nomena including Mach cones, 5 shock waves, 6 tsunamis 
(steepening effect) and solitons. 7 

Monitoring the dust behavior is a dynamic-state es- 
timation problem. Recursive estimation of a dynamic 
state (position/ velocity /acceleration) based on remote 
measurements is often referred to as 'tracking', with the 
tracked object called the 'target'. 8 ' 9 The most widely- 
used technique for tracking individual dusty plasma par- 
ticles based on measurements of their positions is known 
as particle tracking velocimetry (PTV), as used in the 
experiments of Refs. 10-13, for example. In PTV, the 
average velocity of a target over one sample period is de- 
termined by differencing consecutive position measure- 
ments. The instantaneous velocity of a target can be es- 
timated by including a Bayesian inference step in the re- 
cursion. This prediction step is based on a priori knowl- 
edge of the target dynamics. The Kalman filter 8 ' 9 ' 14 is 
an example of a recursive estimator which combines re- 
mote measurements with a predictive step suitable for 
linear dynamics. Nonlinear dynamics can be handled by 
the extended Kalman filter (EKF), which uses a trun- 
cated series expansion of the nonlinear dynamical equa- 
tion (usually first- or second-order). The Kalman filter 



is an optimal estimator for linear systems and Gaussian 
precision, in the sense that the minimum mean-square 
error is obtained. No such result holds for nonlinear sys- 
tems, but the EKF is widely used because it treats the 
nonlinearity explicitly and typically performs very well 
when the initial errors and noises are not too large. 8 A lin- 
ear Kalman filter has been used to track a single particle 
in a simulated dust crystal in Ref. 15. In that work, the 
x- and y-dimensions were filtered independently, which 
limits the accuracy of the prediction. 

In this work we present and implement an EKF-based 
algorithm to track myriad (thousands of) dusty plasma 
particles during the shock wave experiment of Ref. 16, 
and during a computer simulation of a dusty plasma 
shock wave. Central to our algorithm is an interacting 
multiple model 17 ' 18 (IMM) tracker which, using a pre- 
determined set of EKFs, automatically switches to the 
EKF most appropriate to the changeable dust kinemat- 
ics. Synthetic data from the computer simulation (de- 
scribed in Section II) was used to characterize the al- 
gorithm performance in Section III, showing significant 
improvement over PTV. The algorithm was applied to ex- 
perimental data in Section IV, where a pressure-volume 
diagram was generated using a Hugoniot shock-jump re- 
lation. The dusty plasma community is largely unfamil- 
iar with modern target tracking algorithms, that orig- 
inate in aerospace engineering, 8 ' 9 so we have included 
extensive technical details in an appendix. 



II. DUST DYNAMICS 

Each charged dust particle in a dusty plasma creates 
a Coulomb-like potential that is screened by the combi- 
nation of other particles and the ionized gas. The dust 
can be treated as point-sources of charge when the radii 
R are much smaller than the Debye screening length Ad , 
which in turn is smaller than the inter-particle separation 
r. The effective potential experienced by particle j due 



to particle k is a screened, repulsive Coulomb interaction 
of the Debye-Hiickel/ Yukawa form 19 : 

e -r/A D 

0j,fc(fj,fc) = koQl — - — , (1) 

where Qd is the particle charge, fco = l/G^^o) is 
Coulomb's constant, the plasma screening distance is 
the Debye length Ad, and the particles are separated by 
^j,k = r ^j,k- A hat denotes a unit vector. Taking the neg- 
ative of the gradient of the two-particle potential yields 
the effective interaction force between pairs of charged 
dust particles: 

^(^) = ^-~(f- 2 + f- 1 )f J , fc , (2) 

where the dimensionless inter-particle separation is de- 
fined by r = r/Ao- In monolayer experiments, the dust 
is confined by an external potential which typically has 
a shallow parabolic profile (or similar in-plane circular 
geometry). The dust particles align themselves in an im- 
perfect hexagonal lattice. 3,4,7 Such a dust crystal was the 
initial condition for the work in this paper: for both ex- 
periment and simulation. The particles were subjected to 
an external force of short duration that induced a shock 
wave, as in the experiment of Ref. 16 and in the simula- 
tion Run 7 of Ref. 7. The dust was imaged at 1000 frames 
per second in both the experiment and the computer sim- 
ulation (where synthetic images were generated). 

The total force on the jth particle, F, is the sum of 
all two-particle interaction forces F 3 — J2k-Fj,k: plus a 
damping/drag force due to collisions with the plasma, 
plus the force due to the global confinement potential, 
and the short-lived external force (see Ref. 7 for details) . 



Fig. 1. The shock wave traversed the captured field of 
view (the 'scene') in one second, with 1000 simulated 
images generated at millisecond intervals (At = 10 _3 s). 
These images were designed to be faithful replicas of 
those from experiments so that the tracking algorithm 
was tested with realistic data. 
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FIG. 1. Shock wave simulation. Enhanced grayscale image 
of a shock wave traversing a dusty plasma. Bulk/average 
quantities of the dust were calculated for each region of a 
20 x 20 grid. A cross-section 'slice' is highlighted. 



III. ALGORITHM PERFORMANCE 

The performance of any tracking algorithm is char- 
acterized by the errors in the tracks (estimated states). 
A computer simulation is necessary for this because the 
true behaviour of the targets, the 'ground truth', is un- 
known in an experiment. In this section we present errors 
in the dust particle tracks (position and velocity) using 
the ground truth from a molecular dynamics simulation 
of shocked dust. PTV errors were also calculated for 
comparison. Section III A contains an overview of our 
tracking algorithm, with extensive details included as an 
appendix. 

The dynamics of the simulated particles were deter- 
mined using a fifth-order Runge-Kutta numerical inte- 
gration of the Newtonian equations of motion for position 
q(t) and velocity v(t): dq(t)/dt = v(t), dv(t)/dt = F/m, 
where m is the uniform particle mass (see Ref. 7 for pa- 
rameter values used in the simulation). 

We simulated N = 3000 particles in a weakly-confined 
2D dust crystal monolayer, subject to an impulsive ex- 
ternal force that induced a shock wave, as illustrated in 



A. Tracking algorithm 

The general procedure for tracking dust (details in Ap- 
pendix A) starts with particle detection followed by track 
initialization. Each dust particle typically illuminates 
a few pixels in an image. The pixel-intensity- weighted 
centroid of each illuminated region constitutes a particle 
detection, or measurement. For each subsequent image, 
track management (also known as track maintenance) as- 
sociates particle measurements to existing particle tracks, 
or optionally initiates new tracks from unmatched mea- 
surements. Occasionally a particle track will not have 
a corresponding measurement. Such a missed detection 
can occur when a dust particle moves out of the 2D laser 
illumination sheet, or when two particles are close enough 
to illuminate a single region of pixels. After a number of 
consecutive missed detections, a particle track is termi- 
nated. The final stage of the tracking procedure is our fo- 
cus in this paper: state estimation. State estimation uses 
Bayesian inference to predict the target behavior based 
on a physical model describing expected target dynam- 
ics. This prediction is combined with measurement in a 
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weighted sum with weightings assigned so as to minimize 
the expected error in the estimate. For this purpose we 
used a discrete-time extended Kalman filter (EKF), 8 ' 9 
due to the nonlinear nature of the dust interactions (see 
Section II). The first-order EKF is a piecewise lineariza- 
tion of the nonlinear dynamics that is used widely in 
aerospace engineering and works very well in practice. 

Kalman filters (extended or otherwise) recursively up- 
date the state estimate using a weighted average of pre- 
dicted and measured values. The weights are determined 
by the uncertainties in the prediction and measurement 
(stored in covariance matrices). These uncertainties are 
directly influenced by two filter-design parameters known 
as the process noise and the measurement noise, respec- 
tively. Tuning these design parameters allows one to con- 
trollably modify the weighted average between measured 
and predicted values. 

When the dust is excited by an electrostatic force, the 
dynamical model used for prediction is augmented by an 
additional force. In target-tracking nomenclature, this 
temporary deviation from the original dynamical model 
is known as a maneuver. Multiple models can be com- 
bined in an efficient manner using the interacting multi- 
ple model (IMM) approach. 8,17,18 Just as an EKF pro- 
duces a state estimate by combining a predicted value 
with a measured value in a weighted sum, the IMM algo- 
rithm can be thought of as combining multiple state es- 
timates (from multiple EKFs, for example) in a weighted 
sum. The weights for each model in the IMM depend on 
the probabilities of each dynamical model being correct, 
as determined from the measurement-based likelihoods. 
Just as the weights in an EKF are directly affected by 
filter-design parameters (the process noise and measure- 
ment noise), the IMM weights are directly affected by 
selectable parameters known as mode-switching proba- 
bilities. As the name suggests, mode-switching probabil- 
ities affect how quickly the IMM switches between modes 
at the onset and at the conclusion of a maneuver. 

We used an IMM with three prediction models (EKFs) 
for the dust dynamics. The first EKF predicted dust be- 
haviour based on the dust interactions. This is known 
as the base model, or mode, of the IMM. The second 
EKF model added an external force of fixed magnitude 
in the positive-X direction, which accounted for the elec- 
trostatic excitation force that generated the shock. The 
third EKF model extended the base model by adding 
an external force in the negative-X direction, which im- 
proved the tracking of particles reflected off the dust 
cloud bulk. The second and third EKF models are known 
as maneuvering modes. Technical details of the three 
EKFs and the IMM are contained in the appendix in 
sections A3 and A 4, respectively. 

We partitioned the myriad targets into subsets, which 
were tracked independently across multiple computers 
using our EKF-based IMM trackers. This partitioning 
was necessitated by the exponential dependence of com- 
putational load on the number of targets per tracker, AT T : 
the multi-target state vector in two dimensions scales as 



2A^x, with the corresponding covariance matrix scaling 
as the square of this. As discussed in the appendix, this 
would consume (at double precision) in excess of two 
gigabytes of memory per time step to track thousands 
of particles. The subset state estimates were combined 
offline into a single myriad target state, while ensuring 
that particles were counted only once each - there was 
no unchecked overlap of tracks where a single particle 
was accidentally tracked by multiple trackers. Such track 
overlap occur due to particle association errors caused 
by missed detections or intersecting particle trajectories. 
Tracker sizes substantially less than 100 are manageable 
on a desktop computer. Since N T is a design parameter 
in myriad target tracking, we considered iVr = 1, 6 and 
12 particle (s) per tracker for comparison. With the goal 
of minimizing track errors, within limitations imposed by 
the heavy computational resource requirements of multi- 
target tracking, we found that target loss was slightly 
reduced for larger iVr, due primarily to the reduced like- 
lihood of track overlap. Average RMS error levels did 
not vary significantly with TVt- 

B. IMM performance 

The point of the IMM is to detect dust particle ma- 
neuvers. Therefore, we measured the IMM performance 
by comparing the IMM tracker errors with those from a 
single-mode tracker that did not handle maneuvers: an 
EKF tracker based on the Yukawa interactions alone. 

Figure 2 shows a typical result for all tracker configu- 
rations. The bottom plot shows that the IMM trackers 
were a factor of four better than the EKF trackers (1% 
down to 0.25%) at recovering particles lost during the 
confusion created by the shock wave formation between 
t = 0.2s and 0.4s (the maneuver). The top and middle 
plots show that this comes at the cost of increased root- 
mean-square (RMS) errors for estimated position (« 4% 
larger error) and velocity (~ 20% larger error). This tar- 
get loss in the simulation was due to the aforementioned 
track overlap during maneuvers. 

C. Bulk quantity: kinetic energy map 

In Figure 2 we have shown that IMM trackers produce 
significantly higher-accuracy estimates of the dust parti- 
cle kinematics than using PTV, particularly for estimat- 
ing velocity. In this section we use this higher-accuracy 
data to estimate a bulk physical quantity: the average 
kinetic energy 

KE av = . (3) 

Here m is the particle mass, v is the particle speed, and 
the overbar indicates the average over multiple particles 
within each bin in a 20 x 20 grid across the scene (field of 
view). The bins are shown in Fig. 1, which contains the 
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FIG. 2. RMS errors in position (top) and velocity (center), 
and percentage target loss (bottom) as functions of time. 
Tracker configuration: N = 1500 total targets, iV T = 12 PPT. 
Compare the IMM tracker (solid red line) with the base mode 
EKF (dotted blue line) and PTV (dashed green line). 



simulated image at t = 0.499s, enhanced for presentation 
purposes (larger dots and increased brightness). 

A snapshot of the kinetic energy map at t = 0.499s, 
estimated using N — 3000 with TVr = 1, is shown in 
Fig. 3(b), with the corresponding true map in Fig. 3(a) 
showing excellent agreement. Bins containing fewer than 
five particles are ignored and are crossed out in the figure. 
Similar results were achieved with TVt = 6 and 12. 
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FIG. 3. Snapshots of the kinetic energy map at t = 0.499s: 
(a) truth; (b) IMM-tracker estimate using N = 3000 with 
iVr = 1. The outer three rows/columns of empty bins have 
been omitted. Particle locations are overlaid for reference. 
Nonempty bins containing fewer than five particles were ig- 
nored and are crossed out. 

The time-dependence of the accuracy of the estimated 
kinetic energy map can be visualized by combining se- 
lected snapshots of a grid cross-section, or 'slice', taken 
in the direction of shock- wave propagation. We consider 



the central slice highlighted in Fig. 1. Figure 4 contains 
a three-dimensional representation of snapshots of this 
slice taken every 50ms, with the height showing the ki- 
netic energy. Figure 4(a) is the ground truth calculated 
from the simulation, which compares favorably with the 
average kinetic energy estimate for in Figure 4(b). The 
peak of nearly-constant height in the foreground is the 
shock- front. The highest peak (for small X) shows dust 
particles initially excited by the external force, then later 
reflected (at lower speed) off the dust cloud. The kinetic 
energy valley between this excitation/reflection and the 
shock front is where the dust cloud re-crystallizes in the 
wake of the shock wave. 
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FIG. 4. Kinetic energy map cross-section slice (50ms-spaced 
snapshots): (a) Truth (from the simulation); (b) IMM-tracker 
estimate using N = 3000 with Nt = 1. 



IV. EXPERIMENT 

In the experiment reported in Ref. 16, a shock wave 
was generated by a short sequence of negative voltage 
pulses along a wire adjacent to the left side of the dusty 
plasma, and slightly below the plane of the dust cloud. 
The vertical component of this excitation caused many 
dust particles to move out of the laser illumination sheet, 
and hence to disappear, for up to 10 consecutive frames. 
These missed detections were the primary obstacle to 
tracking dust particles in the left of the scene. We tracked 
a total of N = 1860 particles — those that were present 
in all of the first few images. 

Figure 6 shows a snapshot of the experimental results 
at t = 432ms. Shown are the estimated positions and av- 
erage particle number density, velocity, and kinetic tem- 
perature in the direction of the shock wave propagation 
(X). Averages were calculated in each of 50 vertical bins 
across the scene, with the standard deviation for each bin 
shown as error bars (PTV) and a shaded region (IMM 
tracker). The kinetic temperature of particles in each 
vertical bin was calculated as half the particle mass mul- 
tiplied by the mean-squared velocity deviation within the 
bin. Missed detections affected the PTV results severely, 
with many particle velocities (and consequently kinetic 
temperature values) missing — note the large gaps in the 
PTV plots in Fig. 6. Contrast this with the IMM tracker 
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FIG. 5. Side-view schematic of the dusty plasma experi- 
ment. The dusty plasma is contained wholly within a win- 
dowed chamber, with the dust monolayer illuminated by a 
2D laser sheet and imaged from above by a digital camera. 
An offset, below-plane wire (not shown) was used to excite a 
shock wave in the monolayer as in Ref. 6. 



where the problem did not occur. The significant target 
loss which occurred in the shock formation zone (left of 
the scene) due to missed detections was not a problem 
because we were more interested in the shock itself. For 
a recent discussion of problems and limitations of PTV 
using high-speed cameras, we refer the reader to Ref. 20. 

The shock-jump (Hugoniot) relations 21 use conserva- 
tion laws to relate some physical properties of a fluid 
upstream (behind) and downstream (ahead) of an ideal 
shock front. Specifically, the pressure, density and veloc- 
ity are related by three equations that also include the 
shock speed. For known initial conditions downstream, 
and known shock speed, only two of the upstream prop- 
erties need to be determined to find the third unknown 
behind the shock. We consider conservation of momen- 
tum across the shock front, which gives the second shock 
relation 21 



P 2 -P 1 = mn li2 (U s - v^ 2 )(v 2 - v x ) 



(4) 



Here Pi, 2, ^1,2 and v\^ 2 are the pressure, number density 
and particle velocity downstream (1) and upstream (2) 
from the shock. Our IMM tracker results were used to 
determine average bin values for v± j2 and n\. The shock 
speed U a = dX s (t)/dt was determined from the shock 
front position X s (t), which was fitted to the tracker re- 
sults for peak number density nx(i). An ideal shock 
would travel with constant speed c s . Real shocks slow 
down due to damping, which here is a drag force primar- 
ily due to ion/neutral collisions with the dust particles. 
To first order, U s (t) = c s — jt. We determined c S5expt = 
42.7 ± 1.8mm/s, 7 ex pt = H-l ± 1.9mm/s 2 in the experi- 
ment, and c S)S i m = 53.8±1.5mm/s, 7 S i m = 1.7±1.9mm/s 2 
in the simulation. Higher dissipation in the experiment 
7ex P t > 7sim suggests either underestimated dissipation 
in the simulation, or additional unmodelled sources of 
dissipation. Substituting m, ni, vi, v 2 and U a into equa- 
tion (4) yields the shock-induced pressure jump P 2 — P\. 
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FIG. 6. Tracking shocked dust in an experiment: (top to 
bottom) snapshot at t — 432ms of particle positions, particle 
number density, velocity and kinetic temperature. Standard 
deviations for each bin are shown as a grey shadow (IMM) 
and error bars (PTV). Missing error bars are due to missed 
detections where PTV is unable to calculate particle velocity. 



This is plotted in Fig. 7(a) as a function of time, and in 
Fig. 7(b) as a function of n\(t s ) / n 2 {t) — the inverse of 
the shock-induced compression. Here we were restricted 
to times where the shock was present (after formation 
at t s and before leaving the scene): t s < t < 600ms. 
Figure 7(b) is essentially a pressure- volume, or P-V, dia- 
gram. Experiment (crosses, broken line) and simulation 
(triangles, solid line) displayed the same qualitative be- 
havior: shock-induced pressure decreased over time and 
with increasing volume. Pressure decreased over time due 
to the shock wave losing energy as it traversed the dust 
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cloud. Pressure was expected to decrease with increas- 
ing volume, based on basic thermodynamic arguments 
(e.g., Boyle's law for an ideal gas). We speculate that the 
quantitative differences between experiment and simula- 
tion were due in part to dust particle uniformity in the 
simulation versus non-uniformity in the experiment. For 
example, we would expect a non-uniform distribution of 
values for particle mass, charge, etc. in the experiment. 
The large gaps in the PTV velocities shown in Fig. 6 
meant that a P-V diagram could not be calculated using 
PTV. 
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FIG. 7. Shock-induced pressure jump. Experiment (crosses, 
dotted line) and simulation (triangles, solid line) showing a 
steady decrease in shock-induced pressure jump as a function 
of (a) time, and (b) inverse compression. 



data. Estimation of dynamic states based on remote 
measurements is often referred to as target tracking, 8 ' 9 
which here was a nonlinear dynamic estimation problem 
involving myriad maneuvering objects (up to 3000). Ex- 
tensive analysis of the results from computer simulations 
allowed us to quantify the algorithm performance, which 
was significantly more accurate than using particle track- 
ing velocimetry (a factor of three reduction in the aver- 
age error). When applied to the experimental data of 
Ref. 16, the target tracking algorithm demonstrated fur- 
ther superiority to PTV — the tracking algorithm was 
robust to missed detections. Missed detections prevented 
calculation of particle velocities for PTV. Using Hugoniot 
shock-jump relations and tracker-estimated kinematics of 
the dust, the shock-induced pressure jump and compres- 
sion was calculated for the dust cloud. As expected for 
a non-ideal shock, pressure decreased as the shock wave 
lost energy over time. Pressure also decreased as a func- 
tion of increasing volume (inverse compression). Quali- 
tative agreement was observed between experiment and 
simulation. 

In this paper we have shown that the improvement 
in precision and reliability of estimating dust kinematics 
using a target tracking algorithm is important for deter- 
mining the physics of dusty plasmas. Our approach can 
be applied in future experiments to reliably determine 
quantities such as kinetic energy, kinetic temperature, 
diffusion coefficients, 22 and dynamic viscosity. 23 
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Appendix A: State Estimation and Tracking 



V. CONCLUSION 

We have considered a monolayer dusty plasma crys- 
tal disturbed by an electrostatically-induced shock wave, 
presenting results from a computer simulation and the 
experiment of Ref. 16. The dust particle kinematics were 
estimated recursively using an algorithm that combines 
measurements of position (from images) with predicted 
behavior based on multiple models for the particle dy- 
namics. The algorithm provided accurate estimates for 
the motional states of the dust particles (position, ve- 
locity and acceleration) and also selected which of the 
motion models most accurately reflected the measured 



As overviewed in Section III A, estimation of a dynamic 
state consists of combining observations (measurements) 
with predictions. In this appendix we present details 
of the general procedure used in this paper for track- 
ing multiple targets. It consists of the following steps. 
Remote measurements of the targets are made, usually 
at regular intervals. Following track initialization, track 
maintenance occurs after each measurement. Measure- 
ments of the multiple target states (often only positions) 
are either associated to existing tracks, or used to initiate 
new tracks. Tracks lacking a measurement (known as a 
missed detection) are terminated after multiple consecu- 
tive missed detections have occurred. State estimation is 
then implemented and the procedure is repeated for each 
subsequent measurement. 
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1. Measurement: image processing 

Dusty plasma observations consist of a time sequence 
of images taken at (typically) regular intervals separated 
by At. A snapshot at t = 0.353s from the experi- 
ment of Ref. 16 is shown in Figure 8(a). The image 
has been enhanced for display purposes (increased parti- 
cle size and overall image brightness). Similar enhance- 
ment was used on the synthetic image in Figure 8(b). 
Our synthetic images emulate the experimental images 
in Refs. 6, 16, 24, and 25: we also use grayscale, 8-bit, 
tagged image file format (TIFF) images with a 1024-pixel 
square field of view. The field of view, or scene, is 27.3mm 
wide for the experiment, and 80mm wide for the simu- 
lation, corresponding to respective spatial resolutions of 
approximately 26.7/im and 78.1/im per pixel (37.5 and 
12.8 pixels per mm). The wider field of view in the sim- 
ulation enabled the entire dust cloud to be viewed. 

After thresholding the image (setting pixel intensities 
below some threshold equal to zero), a particle is iden- 
tified as a contiguous group of bright pixels. Using the 
moment method, the intensity- weighted centre-of-mass of 
this group of pixels gives the measured/detected position 
of a particle 24,26 : 



Is 



(Al) 



where the intensity of each bright pixel is I 3 , and x 3 is the 
2x1 vector containing the x- and y-pixel number for the 
j th pixel in the group containing N p contiguous bright 
pixels. Here the total intensity of the group of bright 
pixels is = ^2^=1 1 j - For particles illuminating more 
than a few pixels, the moment method yields a typical 
measurement precision of around one tenth of a pixel in 
each direction. More advanced techniques can improve 
this precision, such as those considered in Refs. 24 and 
26, at the cost of heavier computational load. However, 
the potential improvement is limited if measurements are 
not the dominant source of inaccuracies. 




(a) Real image: width 27.3mm, 1024 pixels 




(b) Synthetic image: cropped width « 27.3mm, 351 pixels 



2. Track maintenance 

Maintaining continuous tracks involves a decision pro- 
cess known as 'measurement-to-track association'. Each 
measurement is associated to a track, dismissed as a false 
alarm, or saved as a possible new track. Missed detec- 
tions must also be handled, since object detection can 
never be perfect. For example, dust particles may move 
out of the two-dimensional laser sheet and no longer be 
illuminated. We deleted tracks after nine consecutive 
missed detections. Measurement-to-track association is 
performed online as part of the tracking algorithm. 

Track maintenance is complicated when partitioning 
myriad targets. Unless using centralized, online data fu- 
sion for all trackers (i.e., a parallel tracking algorithm), 



FIG. 8. Enhanced images of a dusty plasma shock wave (a) 
experiment, 16 and (b) simulation. 



track 'overlap' may occasionally occur. This is where a 
single particle is tracked by multiple trackers/computers. 
Track overlap often results when particle trajectories in- 
tersect. Here we delete such duplicate tracks as part of 
our offline data fusion process. 



3. The extended Kalman filter 

We process the measurements and estimate the dust 
particle states (two-dimensional position q, velocity v and 
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acceleration a) with a discrete-time extended Kalman fil- 
ter (EKF). 8 ' 9 The state of each tracked particle at time 
step k = 0,1,2,... is contained in a six-element vector: 
x(k) — (<2x, ^x? <2y? % 5 a y) T • The uncertainty in this 
state is estimated and stored in a corresponding 6 x 6 co- 
variance matrix S$(k). Each step in the filter recursion 
follows measure-update-predict, with the state updated 
after a measurement as in the standard Kalman filter: 

x(k+) = x(k) + K(k) • [z(k) - z(k\k - 1)] . (A2) 

Here the innovation [z(k) — z(k\k — 1)] is weighted by the 
Kalman gain K{k) = S £ (k)-H T -[H-S £ (k)-H T ^R(k)}- 1 . 
The innovation is the difference between the actual mea- 
surement of position z(k) and the expected measurement 
z(k\k — 1) = H • x(k) = (q x (k), q y (k)) T . Here the mea- 
surement matrix is 



H 



1 
1 



(A3) 



The Kalman gain is a matrix version of the ratio of 
track uncertainty to total uncertainty (track plus mea- 
surement), with the measurement uncertainty given by 
ii(fc), the expected covariance matrix for the measure- 
ment. Note that matrix multiplication applies through- 
out — the dots are to guide the eye. 

The updated state x(k+) is now predicted forward to 
the next time step using a discretization of the governing 
dynamical process equations given in Sec. II: 



x(k ■ 



1) = f(x(k+)) 

( 



where the first-order term from the aforementioned Tay- 
lor expansion involves the Jacobian matrix 



q*{k+) + (At)v x (k+) + (At) 2 Fi(x(k + ))/2m 
v x (k+) + (At)FZ(x(k+))/m 
Fi(x(k+))/m 
q y (k + ) + (At)vy(k+) + {AtfF^{x{k + ))/2m 
v y (k+) + {Ai)Fi(x{k + ))/m 
F3(x(k+))/m 

(A4) 



fix) = 



dj_ 
dx 



(A6) 



x=x(kj r ) 



In the standard Kalman filter, the Jacobian is replaced 
by the prediction matrix that represents the linear dy- 
namics. As there is no truncated series expansion in the 
standard Kalman filter, the covariance is a true repre- 
sentation of the prediction error for linear dynamics with 
Gaussian uncertainty. In the EKF, the covariance in the 
prediction is an approximate representation of the pre- 
diction error. The first-order EKF is a piecewise lin- 
earization of the nonlinear particle dynamics. This lin- 
ear approximation should be accurate while the errors 
in the state estimates are small compared to the non- 
linear nature of the underlying dynamical function. 9 For 
our purposes, when the dust particles are in a crystalline 
state, this approximation is excellent. In more dynamic 
moments (during the shock wave, for example), the ap- 
proximation begins to break down and so the tracker per- 
formance is expected to suffer, but not to fail. This is be- 
cause the process noise in Q(k) affords us the ability to 
assign a weight to our confidence in the accuracy of the 
linearized prediction. Similarly, the measurement noise 
in R(k) reflects our confidence in the measurement accu- 
racy. These values can be controllably varied to improve 
the filter accuracy (see Ref. 30 for an example of such 
filter 'tuning'). 

Larger linearization errors can be dealt with by assign- 
ing a lower weight to the prediction, relative to the mea- 
surement. Alternatively, more advanced techniques can 
be used. These include unscented filters 27 and particle 
filters. 28 ' 29 

When tracking multiple targets with a single tracker, 
the single-particle state vectors are stacked to form a 
large vector and the filter matrices are correspondingly 
resized in a block-diagonal fashion. For 3000 particles, 
this would be an 18000-element vector, with correspond- 
ing matrices containing up to 18000 x 18000 elements. At 
double precision, this amounts to approximately 2.5 giga- 
bytes of storage per time step, as well as many floating- 
point operations. For this reason, we partitioned the 
myriad targets into subsets that were manageable as in- 
dividual multi-target trackers. 



where At is the time step. Here we have ignored higher- 
order terms from the Taylor expansion of the full dy- 
namical model equation, which includes the (Gaussian 
random) process noise v(k): x(k + 1) = f(x(k+)) + v(k). 
The process noise assigns a selectable Gaussian uncer- 
tainty to the predicted state /(£), which increases the 
state covariance by Q(fc), the process noise covariance 
(matrix). For the first-order EKF the covariance in the 
predicted state becomes 

S £ (k + 1) = f\x) ■ S A {k+) ■ f(x) T + Q(k), (A5) 



4. Interacting multiple model estimator 

The prediction step inherent in any Kalman-filter- 
based technique for dynamic state estimation is based 
on a physical model for the underlying dynamics of the 
target, as discussed for an EKF in the Section A3. If the 
underlying dynamics change, then the filter performance 
(accuracy) will decrease due to the reduced ability of the 
model to predict the changed/changing dynamics of the 
target. Such dynamical changes, called maneuvers, can 
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often be described by adding an input u(k) to the pre- 
diction step of the filter. Approaches to handling an un- 
known input fall into two broad categories 8 : modeling 
the unknown input as a random process; and estimating 
the unknown input in real time. Input estimation can be 
computationally costly as it involves either maintaining 
a history (sliding window) of estimates that is used to 
detect a maneuver, or it involves increasing the state di- 
mension by augmenting the state vector with the input 
to be estimated. In the case where the underlying system 
dynamics is known to switch between a finite number of 
modes, it can be significantly less costly to run a filter 
for each dynamical model than to perform input estima- 
tion. Such multiple model approaches to state estimation 
use a Bayesian framework where the prior probabilities 
for each model being correct are updated to the poste- 
rior probability at each point in time. In an interacting 
multiple model 17,18 state estimator, the mode probabili- 
ties 11 j are used in conjunction with tunable Markovian 
mode-switching probabilities pij to calculate a mixed 
(weighted) initial condition at each time, which is filtered 
using each possible current model. All of this comes at a 
fraction of the cost of input estimation, with significantly 
improved performance during maneuvers 8 and the added 
benefit of being decision-free (no maneuver detection is 
necessary) . 

The IMM algorithm proceeds as follows. Mixing prob- 
abilities are calculated from the probability Hi(k — 1) 
that mode i was in effect at time step k — 1, given that 
mode j is in effect at time step /c, conditioned on the 
measurements up to k — 1: 

Nj (k-l\k-l) = Pi ^ iik ~ 1] i,j = l...s, (A7) 

C 3 

where Cj ensures normalization and we use s = 3 modes 
(see below). The mode probabilities fij(k) are stored 
along with the state vectors x 3 (k). The mixing probabil- 
ities are used to calculate the mixed initial condition 



(fc-l|fc-l) = ^x i (k-l\k-l)ii i \ j (k-l\k-l), (A8) 



which is input to the j filter, along with the correspond- 
ing covariance 

i=l 

where we have dropped the (k— l\k— 1) arguments for pre- 
sentation purposes. Each filter has a likelihood function 
associated with it, Aj(/c), which is a Gaussian probability 
reflecting the likelihood of obtaining the new measure- 
ment z(k), given the predicted measurement using the 
mixed initial condition, z J (k\k — 1) = H • x 0j ' {k — l\k — 1). 
After running the filters and calculating the likelihood 
functions, the mode probabilities are updated to 



where c is a normalization constant. 

We use a three- mode IMM, with a model each for be- 
fore, during and after the shock wave. The first model 
is an EKF based on our work in Ref. 30, designed for a 
dust crystal. We refer to this as the base mode. The 
second and third modes are referred to as the maneu- 
vering modes, where the target dynamics deviate from 
the base mode (usually for short periods of time). The 
second mode includes an input u(k) to account for an 
external force in the positive-X direction, extending the 
prediction equation (A4) to 



x(k + l) = f(x(k+)) + G y u(k) 

( (At) 2 /2 \ 
At 



/(*(*+)) 



1 



V o / 



i^shockW, (All) 



where G y is the input gain matrix for the unknown ef- 
fective force of tunable magnitude u(k) = -F s hock(fc)- The 
third mode is designed to account for particles reflected in 
the opposite direction, and so the prediction equation has 
the same form as (All), but with u(k) — i 7 'aftershock(^) in 
the negative-X direction (towards the left of the scene). 
The remainder of the EKF equations are unchanged by 
the inclusion of these inputs. For ephemeral maneuvers, 
the base mode should be the best description of the dust 
dynamics for most of the time. Thus, we choose the fol- 
lowing mode-switching probabilities 



0.80 0.10 0.10 
0.30 0.60 0.10 
0.40 0.10 0.50 



(A12) 



Aj(k)cj 



(A10) 



which sum to unity along each row, as required. Larger 
values along the diagonal of pij lead to slower switching 
between modes. This can also be tuned in the opposite 
direction if faster mode switching is desired. Faster mode 
switching tends to yield lower peak RMS error during a 
maneuver, at the expense of greater RMS error during 
quiescent periods. 8 Initially, the mode probabilities are 
equal to one other: fi\ = /12 = ^3 = 1/3. IMM mode 
probabilities can be used to visualize the maneuvers, as 
done for a shocked dusty plasma in Ref. 31. 



5. Range of influence 

The Yukawa force in Equation (2) decreases both expo- 
nentially and polynomially with increasing particle sep- 
aration, thereby reducing the range of influence that 
each dust particle has on other particles. Beyond this 
range, two dust particles exert negligible forces on each 
other. Here we define the range of influence f max rela- 
tive to the nearest neighbor separation tnn as the sep- 
aration at which the two-particle Yukawa force magni- 
tude drops to 1% of the Yukawa force due to the nearest 
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neighbor: i ? J ,/ c (r max )/F J? / c (rNN) = 1%- This occurs for 
f max « 3.81fNN- This range of influence concept was 
used to simplify the particle tracking. The cutoff was 
not used when numerically integrating the equations of 
motion. 
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